//
// Programa que calcula una norma para cada una de las     |
// filas de una matriz bidimensional representada mediante |
// un formato de almacenamiento comprimido.                |
//
// Comando de compilacion: 
//                         gcc -o norma practica2.c -lm
//
#include <assert.h>
#include <limits.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <sys/time.h>

#include "pvm3.h"

#define MAXFILAS 1000 /* Numero de filas maximo que se enviaran del padre a los hijos por cada envio */


void inicializacion(int mytid, double ***pvalores, int **pdimfilas, int metodo, int nproc, int nfilas, int ncols, int *nelementos, int **tids);
void obtener_tids(int *mytid, int *ptid);
void leer_argumentos(int argc, char **argv, int *metodo, int * nproc,int *nfilas, int *ncols);
void paso_parametros(int *tids, int nfilas, int nproc, int mytid, int ptid, int *pfilas, int metodo, int **pos_ini, int **pos_fin, int **dimfilas, double ***valores, int *nelem);
void comprimir_matriz(double ***valores,int *dimfilas,int nfilas,int nproc);
void calcular_normas(double **norma, double **valores, int nfilas,int *dimfilas);
void obtener_norma(double **pnorma, double **valores, int nfilas,int pfilas,int *dimfilas, int mytid, int ptid, int metodo,int nproc, int *tids);
void mostrar_norma(double * norma,int nfilas);
void minmax_vector(int *cadena,int tam);
void guardar_normas(double *norma, int nfilas);

/* muestra un error y finaliza el programa */
void error(char * error) {

	pvm_perror(error);
	pvm_exit();
	exit(0);    
}


/**************************************************************************************************
 * FUNCIONES DE ENVIO/RECIBIR
 **************************************************************************************************/	 


/* envia uno o varios int de un proceso a otro */
void enviar_nint(int *numero, int mytid,int tid, int msgtag, int n) {

	int stride = 1;

	if( pvm_initsend( PvmDataDefault ) < 0 ) 
		error ("ERROR enviar_nint(): pvm_initsend \n");
	if( pvm_pkint(numero, n, stride ) < 0 ) 
			error ("ERROR enviar_nint(): pvm_pkint \n");
	if( pvm_send( tid, msgtag) < 0 ) 
		error ("ERROR enviar_nint(): pvm_send \n");
}


/* 
 * Recibe un int de otro proceso 
 * Se crea un protocolo en el que el primer elemento(int) que se recibe indica cuantos
 * elementos se van a recibir.
 */
void recibir_nint(int *numero[], int mytid, int tid, int msgtag) {

	int stride = 1;	
	int n;                 //numero de elementos a desempaquetar.

	if( pvm_recv( tid, msgtag ) < 0 ) 
		error ("ERROR recibir_nint(): pvm_recv \n");
	if( pvm_upkint(&n, 1, stride ) < 0 ) 
		error ("ERROR recibir_nint(): pvm_upkint \n");		
	if( pvm_upkint(*numero, n, stride ) < 0 ) 
		error ("ERROR recibir_nint(): pvm_upkint \n");
}


/* 
    envia el vector con las dimensiones de las filas que procesara cada proceso
    el primer valor será el numero de valores de dimension que mandamos
 */
void enviar_dimesiones(int tid, int nfilas, int *dimfilas, int salto, int inicio, int fin) {

	int stride = 1;
	int i;

	if( pvm_initsend( PvmDataDefault ) < 0 ) 
		error ("ERROR enviar_dimensiones(): pvm_initsend \n");

	if( pvm_pkint(&nfilas, 1, stride ) < 0 ) 
		error ("ERROR enviar_dimensiones(): pvm_pkint \n");
	
        for(i=inicio; i<=fin; i+=salto){
                if( pvm_pkint(&dimfilas[i], 1, stride ) < 0 ) 
                        error ("ERROR enviar_dimensiones(): pvm_pkint \n");        
        }			

	if( pvm_send( tid, 1) < 0 ) 
		error ("ERROR enviar_dimensiones(): pvm_send \n");
	
}


/*
 *      Los envios se repartiran en grupos de MAXFILAS.
 *      Devuelve el numero de elementos que se mandan al hijo
 */
int enviar_filas(int tid, double **valores, int nfilas, int *dimfilas, int salto, int inicio, int fin ) {

	int stride = 1;
        int nenvios;//Numero de envios de MAXFILAS que se van a efectuar
	int i,j;
	int auxini=inicio;

	int tam = MAXFILAS;// tamano que se va a enviar 
	nenvios = (int)ceil((float)nfilas/(float)MAXFILAS);
	int nelem[nenvios];//Numero de elementos que se van a enviar en cada envio	
	
	int pos=0;
        nelem[0]=0;
        
        int nelementos; // numero de elementos que tendra cada proceso
        
	for(i=inicio; i<=fin; i+=salto){   
	        nelementos += dimfilas[i]; 
                nelem[pos]+=dimfilas[i];
                tam--;
                if(!tam){//Cuando sea la fila 0 no se tiene que aumentar la posicion.
                        tam = MAXFILAS; 
                        pos++;	                    
                        nelem[pos]=0;     
                }
        }tam=MAXFILAS;
   
	for(j=0; j<nenvios; j++){
		if( pvm_initsend( PvmDataDefault ) < 0 )
			error("ERROR enviar_double(): pvm_initsend \n");	
	        if(pvm_pkint(&nelem[j], 1, stride ) < 0 ) 
        		error ("ERROR enviar_filas(): pvm_pkint \n");        			        

                for(i=auxini; i<=fin; i+=salto){

                        if(pvm_pkdouble(valores[i], dimfilas[i], stride ) < 0 ) 
                                error ("ERROR enviar_filas(): pvm_pkdouble \n");	
                                        	    	    
                        free(valores[i]);//Se libera la fila original para que no ocupe memoria.
                        tam--;
                        if(!tam || i==fin){//Cuando sea la fila 0 no se tiene que aumentar la posicion.

                                tam = MAXFILAS;    
                                auxini = i+salto;
                                if( pvm_send( tid, j+3) < 0 ) 
                                    error ("ERROR enviar_dimensiones(): pvm_send \n"); 
                                break;   
                        }                        
                }
        }	
        
        return nelementos;
}


// recibe las filas por orden segun su mstag
void recibir_filas(double ***valores, int *pfilas, int *dimfilas, int mytid, int ptid) {

	int aux = 0;
	int nenvios;  //numero de envios que se recibiran
	int i;
	int tam;      //Donde se desempaquetara el tamano del envio
	int stride =1;

	nenvios = (int)ceil((float)*pfilas/(float)MAXFILAS);	
	
	for(i=0;i<(*pfilas);i++)
		aux+=dimfilas[i];		

	//Creacion vector en el que se guardara la matriz
	*valores =(double **) malloc((*pfilas) * sizeof(double *));
	(*valores)[0]=(double *)malloc(aux * sizeof(double));

	aux=0;
	for(i=0; i<nenvios; i++){
		if( pvm_recv(ptid, 3+i ) < 0 ) 
			error ("ERROR recibir_filas(): pvm_recv \n");
		if( pvm_upkint(&tam, 1, stride ) < 0 ) 
			error ("ERROR recibir_filas(): pvm_upkint \n");	
		if( pvm_upkdouble((valores[0][0])+aux, tam, stride ) < 0 ) 
			error ("ERROR recibir_filas(): pvm_upkint \n");		
		aux+=tam;
	}
	
	aux= dimfilas[0];
	for(i=1;i<(*pfilas);i++){
		(*valores)[i] = (valores[0][0])+aux;
		aux+=dimfilas[i];
	}		        
}

// empaqueta el vector de normas y como primer componente manda el tamaño
void enviar_normas(double *norma,int ptid,int pfilas) {
		
		int stride = 1;

		if( pvm_initsend( PvmDataDefault ) < 0 )
			error("ERROR enviar_normas(): pvm_initsend \n");	
		if(pvm_pkint(&pfilas, 1, stride ) < 0 ) 
        		error ("ERROR enviar_normas(): pvm_pkint \n");         
	        if(pvm_pkdouble(norma, pfilas, stride ) < 0 ) 
        		error ("ERROR enviar_normas(): pvm_pkdouble \n");          		
		if( pvm_send( ptid, 2) < 0 ) 
			error ("ERROR enviar_dimensiones(): pvm_send \n"); 			
}


// recibe las normas y las desempaqueta en la posicion correcta del vector
void recibir_normas(int *tids,double **pnorma,int nfilas,int nproc,int metodo,int pfilas) {
	
	int i;
	int stride = 1;
	int nnorma;        // tamaño del vector norma
	int aux = pfilas;
	
	for (i = 1; i < nproc; i++) {

		if( pvm_recv(tids[i], 2 ) < 0 ) 
			error ("ERROR recibir_filas(): pvm_recv \n");		
		if( pvm_upkint(&nnorma, 1, stride ) < 0 )  
			error ("ERROR recibir_filas(): pvm_upkint \n");	
		if (!metodo) {  //consecutivo
			if( pvm_upkdouble((*pnorma)+aux, nnorma, 1 ) < 0 ) 
				error ("ERROR recibir_filas(): pvm_upkint \n");	
			aux += nnorma;		
		}		
		else 		//ciclico
			if( pvm_upkdouble((*pnorma)+i, nnorma, nproc ) < 0 ) 
				error ("ERROR recibir_filas(): pvm_upkint \n");		

	}
	
}

/*************************************************************************************************/
/* MAIN                                                                                          */
/*************************************************************************************************/
int main( int argc, char *argv[]) {

	double **valores, *norma;
	int *dimfilas;
	int nfilas, ncols, nelementos, metodo, nproc;
	int mytid, ptid;
	int *tids;	

	int pfilas; /* filas que analizara cada proceso */	
	int *pos_ini, *pos_fin;
	

	
        struct timeval t0, t1, t;
	
	obtener_tids(&mytid, &ptid);
			
	if(ptid < 0) {
		leer_argumentos(argc, argv, &metodo, &nproc, &nfilas, &ncols);
                inicializacion(mytid, &valores, &dimfilas, metodo, nproc, nfilas, ncols, &nelementos, &tids);	
        }
        
	int *nelem;
	nelem = (int *)malloc(sizeof(int) * nproc);
	nelem[0] = nelementos;        

        //obtener tiempo inicial
        assert (gettimeofday (&t0, NULL) == 0);

	paso_parametros(tids, nproc, nfilas, ptid, mytid, &pfilas, metodo, &pos_ini, &pos_fin, &dimfilas,&valores,nelem);
	
	if ((ptid < 0) && (metodo))  // comprimimos la matriz para que en el padre las filas esten consecutivas en el ciclico
                comprimir_matriz(&valores, dimfilas, nfilas, nproc);  

	obtener_norma(&norma, valores, nfilas, pfilas, dimfilas, mytid, ptid, metodo, nproc, tids);
	
	// obtener tiempo final
	assert (gettimeofday (&t1, NULL) == 0);
	timersub(&t1, &t0, &t);
	
	if (ptid < 0){
		//mostrar_norma(norma,nfilas);
		printf ("Tiempo en paralelo     = %ld:%ld(seg:mseg)\n", t.tv_sec, t.tv_usec/1000);
		// como el primer elemento tiene el total le restamos los que se llevan los hijos 
		int i;
		for (i = 1; i < nproc; i++) 
			nelem[0] -= nelem[i];

		minmax_vector(nelem,nproc);		
		guardar_normas(norma, nfilas);
	}
	
	/* al salir del programa salimos de pvm */
	if( pvm_exit() < 0 ) {
		pvm_perror ("ERROR main(): pvm_exit \n");
	}
	return 0;
}

/*************************************************************************************************/

// funcion que compara dos enteros y devuelve el menor
int min(int a, int b) {

        if(a<b) 
                return a;
        return b;
}

// guarda el vector de  normas en un fichero
void guardar_normas(double *norma, int nfilas) {
        FILE *f; 
        int i;
        
        f = fopen("./normas_paralelo.txt","w");
        for(i=0; i<nfilas; i++){
                fprintf(f,"%f\n",norma[i]);
        }
}

// calcula el mayor y menor elemento de una cadena (cadena) de tamaño (tam)
void minmax_vector(int *cadena,int tam) {

	int i;
	int min,max;

	min = cadena[0];
	max = cadena[0];

	for (i = 1; i < tam; i++) {
		if (cadena[i] > max) max = cadena[i];
		if (cadena[i] < min) min = cadena[i];
	}     
	
	printf("Elementos minimos [%d]\n",min);
	printf("Elementos maximos [%d]\n",max);
	printf("Diferencia      [%d]\n",max-min);

}


// obtiene el id actual y el del proceso padre en pvm 
void obtener_tids(int *mytid, int *ptid) {

	(*mytid) = pvm_mytid();	
	(*ptid)  = pvm_parent();	
}


// muetra los valores del vector norma
void mostrar_norma(double * norma,int nfilas) {

	int i;

	printf("NORMA\n");
	for (i = 0; i < nfilas; i++) {
		printf("%f\n",norma[i]);
	}
}


// obtiene los argumentos del programa
void leer_argumentos(int argc, char **argv, int *metodo, int *nproc, int *nfilas, int *ncols) {	

	if (argc != 5) {
		printf ("Uso: metodo nproc nfilas ncols \n");
		exit(0);
	}
	*metodo = atoi(argv[1]);
	*nproc  = atoi(argv[2]);
	*nfilas = atoi(argv[3]);
	*ncols  = atoi(argv[4]);
}

// en el caso ciclico mete en las primeras posiciones las filas procesadas
void comprimir_matriz(double ***valores,int *dimfilas,int nfilas,int nproc) {
        
        int i;
        int j = 1;
        
        for (i = nproc; i < nfilas; i += nproc) {
                (*valores)[j] = (*valores)[i];
                dimfilas[j] = dimfilas[i];
                j++;
        }
}


// en el caso ciclico cambia de posicion los valores de la norma segun el salto
void descomprimir_vector(double *valores,int nfilas,int salto) {
        
        int i;
        
        for (i = nfilas-1; i >= 0; i--) 
		valores[salto*i] = valores[i];
}


// Inicializa la matriz dispersa
void crear_matriz(double ***pvalores, int **pdimfilas, int metodo,int nfilas, int ncols, int *nelementos) {

	int ncolsfila;
	int i,j;

	if( ((*pdimfilas)=(int *) malloc(sizeof(int)*nfilas)) == NULL ) 
		error ("error: memoria insuficiente para dimfilas\n");

	if( ((*pvalores)=(double **) malloc(sizeof(double)*nfilas)) == NULL )//sizeof(double *) ?
		error ("error: memoria insuficiente para valores\n");

	for (i=0; i<nfilas; i++) {
		ncolsfila = (int)(rand()%(9*i+1));
		if ( ncolsfila > ncols ) ncolsfila=ncols;
		(*pdimfilas)[i] = ncolsfila;
		if( ((*pvalores)[i]=(double *) malloc(sizeof(double)*(*pdimfilas)[i])) == NULL ) 
			error("error: memoria insuficiente para valores\n");
	}

	(*nelementos) = 0;

	for (i=0; i<nfilas; i++) {
		for (j=0; j<(*pdimfilas)[i]; j++) {
			(*pvalores)[i][j] = (double) ( (double) random() / (double) RAND_MAX );
		}
		(*nelementos) += (*pdimfilas)[i];
	}	
}

// crea los procesos hijo
void crear_hijos(int mytid, int nproc, int **tids) {

	(*tids)[0] = mytid;		

	/* crea procesos hijo */
	if( (pvm_spawn("p2", (char**)0, 0, "", nproc-1, &(*tids)[1])) < nproc-1 ) 
		error ("ERROR crear_hijos(): Not enough processes created \n");

}


/*
 *    Inicializacion de la matriz y creacion de los procesos hijos
 */
void inicializacion(int mytid, double ***pvalores, int **pdimfilas, int metodo, int nproc,int nfilas, int ncols, int *nelementos, int **tids) {

        // inicializar vector de procesos
        if ( ((*tids)=(int *)malloc(nproc*sizeof(int))) == NULL ) 
	        error ("ERROR inicializacion(): Not enough memory.\n");	

        crear_matriz(pvalores, pdimfilas, metodo, nfilas, ncols, nelementos);
        crear_hijos(mytid,nproc,tids);
		
}


// calcula los vectores L y U
void calcular_posiciones(int **pos_ini,int **pos_fin, int metodo,int nfilas,int nproc, int *salto) {
	int i;	
	int np;

	if((*pos_ini = (int *)malloc(nproc*sizeof(int))) == NULL)
		error ("ERROR calcular_posiciones(): Not enough memory.\n");	
	if((*pos_fin = (int *)malloc(nproc*sizeof(int))) == NULL)
		error ("ERROR calcular_posiciones(): Not enough memory.\n");		

	if(!metodo){  /* consecutivo */					
                *salto = 1;		
		np =(int) ceil((float)nfilas/(float)nproc);	
		for(i=0;i<nproc;i++){		                        	        
			(*pos_ini)[i]= np*i;
		        (*pos_fin)[i]= min((i+1)*np-1, nfilas-1);
		}			
	}else{        /* ciclico */
	        *salto = nproc;
		for(i=0;i<nproc;i++){		                        	        
			(*pos_ini)[i]= i;
		        (*pos_fin)[i]= nfilas-1;
		}				
	}
}


// funcion que se encarga de enviar los distintos parametros necesarios para los procesos hijo
void paso_parametros(int *tids, int nproc, int nfilas, int ptid, int mytid, int *pfilas, int metodo, int **pos_ini, int **pos_fin, int **dimfilas, double ***valores, int *nelem) {
		
	int i;
	int salto;
	int nfilas_p[2];
	int fin;
	
	if(ptid < 0){ //padre
		calcular_posiciones(pos_ini, pos_fin, metodo, nfilas, nproc, &salto);//rellena salto
		/* Envio numero filas */
        	*pfilas = (int)floor((float)((*pos_fin)[0]-(*pos_ini)[0])/salto) +1;//Numero de normas a realizar por el padre
		for(i=1; i<nproc; i++){		                			        
		        nfilas_p[0] = 1;
                        /* numero de filas a cada proceso floor((pos_fin[i]-pos_ini[i])/salto)+1 */	        
		        nfilas_p[1] =  (int)floor((float)((*pos_fin)[i]-(*pos_ini)[i]) /salto) +1;
                        enviar_nint(nfilas_p, tids[0], tids[i], 0, 2);
                        if(nfilas_p[1] < 1)
                        	break;//A partir de aqui todos los siguientes envios serian de 0 filas
                        enviar_dimesiones( tids[i], nfilas_p[1], *dimfilas, salto, (*pos_ini)[i], (*pos_fin)[i]);
                        
                        // calculo del numero de la parada para enviar filas
                        if (!metodo) fin = (*pos_fin)[i];
                        else fin = (*pos_ini)[i]+(nfilas_p[1]-1)*nproc;

                        nelem[i] = enviar_filas(tids[i], *valores, nfilas_p[1],*dimfilas, salto, (*pos_ini)[i], fin);
		}
				
	}else{	//hijo
		recibir_nint(&pfilas, mytid, ptid, 0);
		if(*pfilas < 1)
			return;
	     
		*dimfilas = malloc(sizeof(int) * (*pfilas)); // vector de dimensiones
		recibir_nint(dimfilas, mytid, ptid, 1);

		recibir_filas(valores, pfilas, *dimfilas, mytid, ptid);		
	}
}


/*
 * Funcion que calcula las normas de la matriz valores dada
 * Este fragmento de codigo sera ejecutado por todos los procesos del programa paralelo
*/
void calcular_normas(double **pnorma, double **valores, int pfilas,int *dimfilas) {
 
        int i, j;
        double normafila;

	for (i=0; i<pfilas; i++) {
		normafila = 0;
		for (j=0; j<dimfilas[i]; j++) {
			normafila += valores[i][j] * valores[i][j];
		}
		(*pnorma)[i] = sqrt(normafila);
	}
	
}

// hace el calculo de las normas y recompone el vector de normas final
void obtener_norma(double **pnorma, double **valores, int nfilas,int pfilas,int *dimfilas, int mytid, int ptid, int metodo,int nproc, int *tids) {
	
	if( ((*pnorma)=(double *) malloc(sizeof(double)*pfilas)) == NULL ) {
		printf ("error: memoria insuficiente para norma\n");
		exit(-1);
	}

	calcular_normas(pnorma, valores , pfilas, dimfilas);
	
	if (ptid<0) { //padre
	        // obtenemos memoria para todas las normas de todos los hijos
		if( ((*pnorma)=(double *)realloc((*pnorma),sizeof(double)*nfilas)) == NULL ) {
			printf ("error: memoria insuficiente para norma\n");
			exit(-1);
		}
		if (metodo) // en el caso ciclico
			descomprimir_vector(*pnorma,pfilas,nproc);

		recibir_normas(tids,pnorma,nfilas,nproc,metodo,pfilas);
				
	 }else //hijo
		enviar_normas(*pnorma,ptid,pfilas);
}


